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ABSTRACT 

A  theoretical  application  of  transient  queueing  analysis  is  provided  for  military  air  traffic  control.  The 
exact  distribution  of  the  nth  arriving  or  departing  flight’s  sojourn  time  in  an  M/M /s  queue  with  k 
flights  initially  present  is  reviewed.  Algorithms  previously  developed  for  computing  the  covariance 
between  sojourn  times  for  an  M/M / 1  queue  with  k  >  0  flights  present  at  time  zero  are  provided  and 
utilized.  Maple  computer  code  is  utilized  for  practical  applications  in  air  traffic  control  of  transient 
queue  analysis  for  many  system  measures  of  performance  without  regard  to  traffic  intensity  (i.e.,  the 
system  may  be  unstable  with  traffic  intensity  greater  than  one),  thus  negating  the  need  for  simulation. 

1  INTRODUCTION 

Many  traditional  simulation  studies  analyze  queueing  systems  in  steady-state,  requiring  appropriate 
warm-up  periods  and  associated  long  simulation  runs.  However,  in  many  cases  the  system  being 
modeled  never  reaches  steady-state;  thus  steady-state  simulation  results  do  not  accurately  portray  the 
system  behavior,  as  is  often  the  case  in  military  air  traffic  control.  The  ability  to  analyze  transient 
results  associated  with  such  models  is  often  complicated  by  intractable  theory,  leaving  simulation 
as  the  only  method  for  analysis.  Further  complicating  the  transient  analysis  is  the  effect  of  initial 
conditions  (Kelton  and  Law,  1985).  Since  steady-state  results  depend  on  running  the  system  long 
enough  to  mitigate  the  impact  of  initial  conditions,  these  steady-state  results  reveal  nothing  about  the 
transient  behavior  of  the  queueing  system.  Our  purpose  here  is  to  combine  new  and  existing  results 
in  transient  queueing  analysis  with  a  symbolic  engine  in  computational  probability. 

There  are  many  classes  of  queueing  systems  where  a  transient  analysis  is  required,  e.g.,  military 
applications  often  use  queueing  models  that  never  reach  equilibrium.  Recognizing  the  need  to  develop 
theory  for  transient  results,  as  opposed  to  steady-state  results,  has  resulted  in  a  wide  literature  in  this 
area.  Initial  work  in  transient  analysis  ironically  appeared  as  an  attempt  to  measure  when  a  system 
achieved  equilibrium.  Law  (1975)  notes  the  consequences  of  failing  to  adequately  account  for  the 
initial  transient  period,  leading  to  Gafarian  et  al.  (1976)  outlining  a  comprehensive  framework  for  the 
initial  transient  problem.  Morisaku  (1976)  addresses  the  time  to  equilibrium  in  simulations  modeling 
the  M/M / 1  queue  and  provides  schematics  for  the  transition  probabilities  given  k  >  0  customers 
initially  present  at  time  t  =  0.  Grassmann  (1977)  compares  three  methods  for  finding  transient 
solutions  in  Markovian  queueing  systems,  Runge-Kutta,  Liou’s  method,  and  randomization,  where 
randomization  is  shown  as  superior  for  large  sparse  transition  matrices.  Pegden  and  Rosenshine 
(1982)  provide  a  closed-form  solution  for  the  probability  of  exactly  i  arrivals  and  j  servicings  over 
a  time  horizon  of  length  t  in  an  M/M / 1  queue  starting  empty  and  idle,  allowing  the  calculation  of 
certain  performance  measures  for  a  specified  time  period.  Odoni  and  Roth  (1983)  take  an  empirical 
approach  to  compare  observed  and  predicted  transient  state  queue  length  for  the  M/M/ 1  queue,  noting 
that  for  small  values  of  t  the  expected  queue  length  is  strongly  influenced  by  initial  conditions,  and 
provide  a  good  approximation  for  an  upper  bound  of  time  to  steady-state.  Kelton  and  Law  (1985) 
consider  the  M/M /s  (s  >  1)  queue  and  provide  expressions  for  calculating  the  probabilities  of  having 
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up  to  n  +  k  customers  in  the  system  upon  the  arrival  of  the  nth  customer,  where  k  is  the  number  of 
customers  in  the  system  at  time  t  —  0.  They  then  apply  these  calculations  to  a  variety  of  measures  of 
performance  with  implications  to  convergence  on  steady-state  delays  and  offer  methods  for  choosing 
queue  initialization  in  simulation.  Much  of  the  work  in  this  paper  is  motivated  by  their  results.  Kelton 
(1985)  extends  the  previous  work  by  considering  M/Em/ 1  and  Em/M / 1  queues.  Parthasarathy  (1987) 
provides  a  transient  solution  for  the  probability  that  there  are  n  customers  in  the  system  at  time  t  for 
an  M/M / 1  queue.  Abate  and  Whitt  (1988)  use  Laplace  transformations  to  analyze  some  transient 
results  of  interest  in  the  M/M/ 1  queue.  Leguesdron,  et  al.  (1993)  provide  transient  probabilities  for 
the  M/M / 1  queue  by  inverting  the  generating  function  of  the  uniformized  Markov  chain  describing 
the  M/M / 1  process.  Transient  distributions  of  cumulative  reward  are  calculated  by  e  Silva,  et  al. 
(1995)  using  uniformization,  where  the  distribution  of  cumulative  reward  is  over  a  finite  interval  with 
reward  rates  represented  by  Markov  model  states.  Grassmann  (2008)  investigates  warm-up  periods 
in  simulation  and  shows  that  in  many  cases  these  warm-up  periods  should  not  be  used,  especially  if 
the  simulation  begins  in  a  high  probability  state. 

The  purpose  of  this  paper  is  to  adapt  previous  work  on  transient  analysis  to  air  traffic  control. 
Kaczynski  et  al.  (2010)  provide  an  extensive  literature  review  in  transient  analysis.  In  this  paper  we 
will  focus  on  the  transient  analysis  of  the  M/M / 1  and  the  more  general  M/M /s  queues  as  related 
to  air  traffic  control,  specifically  on  the  distribution  of  the  nth  flight’s  (arriving  or  departing)  sojourn 
time,  which  is  the  sum  of  the  nth  flight’s  delay  time  (referred  to  as  holding  time  in  aviation)  and  service 
time.  Almost  all  the  references  listed  in  Kaczynski  et  al.  (2010)  address  measures  of  performance 
specified  over  a  finite  time  interval  and  are  the  results  of  numerical  work,  which  is  in  stark  contrast 
to  the  measures  proposed  here,  which  are  based  on  the  exact  distribution  of  specific  flights.  Rather 
than  arriving  at  the  results  numerically,  a  computer  algebra  system  calculates  exact  measures  of 
performance  based  on  a  given  number  of  arriving  and  departing  flights. 

The  M/M/s  queue  is  defined  in  Section  2  for  a  positive  integer  s,  and  a  method  is  given  for 
calculating  the  probability  distribution  of  the  number  of  aircraft  an  arriving  or  departing  flight  sees 
upon  contact  with  air  traffic  control  (ATC)  (upon  arrival  to  an  M/M/s  queue).  Section  3  describes 
how  the  sojourn  time  distribution  is  calculated  for  a  given  flight  in  an  M/M/s  queue  with  k  flights 
initially  present  in  the  system,  k>0.  Section  4  includes  examples  using  the  implemented  procedures  to 
calculate  exact  sojourn  time  distributions,  related  measures  of  performance,  and  graphical  illustrations 
for  varying  parameters  such  as  traffic  intensity  and  number  of  flights  in  the  system.  Section  5  applies 
two  approaches  for  calculating  the  covariance  and  correlation  among  flights  in  an  M/M/ 1  queue. 
Section  6  extends  the  covariance  and  correlation  calculations  by  automating  the  process  of  finding 
the  joint  probability  distribution  function  between  two  flights,  and  provides  the  exact  covariance  and 
correlation  calculations  for  varying  traffic  intensities.  Section  7  concludes  the  paper  by  reviewing 
the  content  and  offering  areas  of  further  study.  Commented  code  is  available  for  all  computations 
conducted  here  upon  request. 

2  BASICS  OF  THE  M/M/s  QUEUE 

The  M/M/s  queue  is  governed  by  iid  exponential  interarrival  times  (the  arrival  stream  is  a  Poisson 
process)  with  arrival  rate  ,  and  iid  exponential  service  times  among  s  identical  servers,  each  with 
service  rate  .  The  interarrival  times  and  the  service  times  are  mutually  independent.  The  traffic 
intensity  of  the  system  is  =  js  .  The  system  consists  of  a  single  queue  with  customers  waiting 
to  be  serviced  by  one  of  the  identical  s  parallel  servers.  If  an  arriving  customer  finds  at  least  one  idle 
server,  the  customer  immediately  proceeds  to  service;  otherwise  the  customer  joins  the  single  queue 
of  those  waiting  for  service  in  a  first-come,  first-served  manner.  To  achieve  classic  steady-state  results 
the  traffic  intensity  must  satisfy  <  1.  This  critical  assumption  is  not  required  in  transient  analysis, 
described  here,  because  the  system  of  interest  never  reaches  equilibrium.  In  this  application,  if  .v  =  1, 
we  will  assume  there  is  a  single  controller  and  a  single  runway.  If  s  >  1,  we  will  assume  there  is 
one  runway  for  each  controller  and  that  the  controllers  have  adequate  space  to  independently  serve 
flights. 

Let  /\(n,f )  be  the  probability  that  upon  the  arrival  of  the  nth  flight  there  are  exactly  i  flights 
in  the  system  including  the  nth  flight  (in  queue  or  in  service),  given  k  flights  are  present  at  time 
t  =  0.  Using  propositions  provided  by  Kelton  and  Law  (1985)  and  a  recursion  algorithm,  [\(n.  i)  for 
i  =  1,2, . . .  ,n  +  k  can  be  computed.  Using  these  probabilities,  it  is  possible  to  find  the  distribution  of 
the  sojourn  time  for  the  nth  flight  in  an  M/M/s  queue,  given  k  flights  are  present  at  time  t  =  0. 

Kaczynski etal.(2010)includestheseresultscodedintheMapleprocedureQueue  (X,  Y,  n,  k,  s), 
where 
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•  X  is  the  exponential  interarrival  time  distribution, 

•  Y  is  the  exponential  service  time  distribution, 

•  n  is  the  index  of  the  flight  of  interest, 

•  k  is  the  number  of  flights  in  the  system  at  time  t  =  0, 

•  s  is  the  number  of  identical  parallel  servers. 

The  procedure  is  written  in  Maple  and  uses  A  Probability  Programming  Language  (APPL),  which 
can  be  downloaded  for  free  at  www .  APPLsoftware  .  com  and  is  described  in  Glen  et  al.  (2001). 
We  choose  to  calculate  the  distribution  of  the  sojourn  time  because  it  is  a  purely  continuous  random 
variable  enabling  us  to  exploit  associated  procedures  in  APPL. 

3  CREATING  THE  SOJOURN  TIME  DISTRIBUTION 

Once  the  necessary  Pk(n, i),  i  —  1,2 ,...,n  +  k,  probabilities  are  calculated,  the  exact  sojourn  time 
distribution  for  the  nth  flight  can  be  calculated.  We  define  Xn  as  the  number  of  flights,  including 
flight  n,  in  the  system  at  time  t,  the  arrival  time  of  the  nth  flight.  The  possible  values  of  X,,  can  vary 
from  a  minimum  of  1 ,  which  occurs  when  flight  n  arrives  to  empty  airspace,  to  a  maximum  of  n  +  k, 
which  occurs  when  0  completions  occur  prior  to  flight  n’s  arrival.  The  mathematical  derivations  for 
both  the  M/M/ 1  and  M /M/s  queues  make  extensive  use  of  the  memoryless  property,  permitting  the 
construction  of  the  distribution  of  Tn,  the  sojourn  time  of  flight  n.  We  present  each  case  separately 
below. 

3.1  Distribution  of  Tn  for  the  M/M/1  queue 

Kaczynski  et  al.  (2010)  show  the  M/M/ 1  queue  sojourn  time  distribution  for  k  =  0  initial  flights 
generalizes  very  elegantly  to  include  k  >  0,  as  indicated  in  Table  1 .  Line  i  of  the  table  occurs  with 
probability  Pk[n,  i)  and  lists  the  distribution  of  the  sojourn  time  for  the  nth  flight,  conditioned  on  i 
flights  being  in  the  system  upon  its  arrival. 

Table  1:  Conditional  sojourn  time  distributions  for  the  M/M/1  queue. 


X„ 

Delay 

Service 

Conditional  sojourn  time  distribution 

1 

0 

exponential  ) 

exponential  ) 

2 

exponential  ) 

exponential  ) 

Erlang!  ,  2) 

3 

Erlangt  ,  2) 

exponential  ) 

Erlang!  ,  3) 

4 

Erlang(  ,  3) 

exponential  ) 

Erlang!  ,  4) 

n  +  k 

Erlang(  ,  n  +  k—  1) 

exponential  ) 

Erlang!  ,  n  +  k) 

Let  gi(t)  be  the  PDF  of  an  Erlang(  ,  i)  random  variable.  Using  the  conditional  sojourn  time 
distributions  for  i  =  1,2 ,...,n  +  k  potential  flights  in  the  system,  each  with  probability  Pk(n,i),  the 
probability  density  function  (PDF)  for  the  nth  flight’s  sojourn  time  T„  is  the  mixture 

n+k 

fn{t)=  Pk(n,i)gi(t)  t>  0.  (1) 

i—  1 


3.2  Distribution  of  Tn  for  the  M/M/s  queue 

Given  s  >  1  parallel  identical  servers  (controllers),  the  nth  flight’s  sojourn  time  distribution  is  still  a 
mixture  of  n  +  k  conditional  sojourn  time  distributions.  However,  each  distribution  might  be  more 
complicated.  For  illustration,  consider  an  M/M/3  queue  starting  empty  and  idle  with  exponential!  ) 
arrivals  and  three  identical  exponential!  )  servers.  It  is  clear  that  for  flights  1,  2,  and  3,  the  sojourn 
time  is  exponential  )  since  all  three  flights  proceed  directly  to  service.  Therefore,  in  the  general 
case,  for  the  number  of  flights  in  the  system  including  flight  n,  which  we  defined  as  Xn,  when  Xn  <  s 
the  conditional  sojourn  time  distribution  is  exponential!  )•  However,  if  Xn  >  s,  then  the  nth  flight 
experiences  a  delay  while  observing  Xn  —  s  service  completions.  When  s  >  1  and  Xn  >  s ,  the  service 
distribution  observed  by  flights  in  queue  is  exponential  with  rate  s  .  Using  this  result,  it  is  apparent 
that  the  delay  time  for  the  nth  flight  is  the  sum  of  A,,  —  s  independent  exponential^  )  random  variables, 


1397 


Kaczynski,  Leemis  and  Drew 


and  using  (1)  is  Erlang©  ,X„  — s ).  To  calculate  the  nth  flight’s  sojourn  time  for  a  particular  value 
of  Xn,  we  sum  the  delay  time  and  the  service  time.  Table  2  shows  the  distributions  conditioned  on 
the  number  of  flights  Xn  encountered  by  flight  n  (including  itself)  for  the  M/M/3  queue,  given  k  =  0 
flights  present  at  time  t  =  0.  The  APPL  procedure  Convolution  calculates  the  distribution  of  a 
sum  of  independent  random  variables.  We  use  the  symbol  ©  to  represent  convolution. 

Table  2:  Conditional  sojourn  time  distributions  for  the  M/M/3  queue  with  k  =  0. 


Xn 

Delay 

Service 

Conditional  sojourn  time  distribution 

1 

0 

exponential  ) 

exponential!  ) 

2 

0 

exponential  ) 

exponential!  ) 

3 

0 

exponential  ) 

exponential!  ) 

4 

exponential(3  ) 

exponential  ) 

exponential(3  )  ©  exponential!  ) 

5 

Erlang(3  ,  2) 

exponential  ) 

Erlang(3  ,  2)  ©  exponential!  ) 

n 

Erlang(3  ,  n  -  3) 

exponential  ) 

Erlang(3  ,  n  —  3)  ©  exponential!  ) 

Since  Xn  represents  the  number  of  flights  in  the  system  upon  arrival  of  the  /7th  flight,  including 
itself,  the  first  row  in  Table  2  corresponds  to  flight  n  arriving  to  an  empty  system  and  the  last  row 
corresponds  to  no  service  completions  prior  to  flight  n’s  arrival.  The  general  form  for  the  M /M/s 
sojourn  time  PDF  is  identical  to  (1),  however,  in  the  M /M/s  case  each  gi(t)  can  potentially  require 
an  additional  step  to  calculate  the  distribution  of  a  sum  of  random  variables. 

4  AIR  TRAFFIC  CONTROL  TRANSIENT  ANALYSIS  APPLICATIONS 

It  is  apparent  that  calculating  the  conditional  sojourn  time  distribution  for  large  n  is  tedious.  Kelton 
and  Law  (1985)  acknowledge  the  computational  difficulty  in  achieving  the  Pk(n,i )  probabilities 
alone.  Conducting  the  added  steps  of  up  to  n  —  s  convolutions  for  the  M/M/s  queue  and  then 
mixing  the  resulting  conditional  distributions  with  the  appropriate  probabilities  can  be  complicated  to 
implement.  Due  to  these  computational  difficulties,  simulation  is  often  required  for  analysis,  however, 
APPL  provides  the  underlying  computational  engine  to  achieve  exact  results  for  such  problems.  As 
mentioned  earlier,  the  APPL  procedure  Queue  (X,  Y,  n,  k,  s)  returns  the  exact  sojourn  time 
distribution  for  flight  n.  APPL  is  capable  of  symbolic  results,  as  illustrated  in  Examples  1  and  2. 

Example  1.  Consider  an  M/M/ 1  queue  (a  single  controller,  single  runway  facility) 
with  arrival  rate  and  service  rate  starting  empty  and  idle  at  time  t  =  0.  For  the 
fourth  flight,  calculate  the  probabilities  Pq(4.  i)  for  i  =  1,2, 3, 4. 

The  APPL  command  MMsQprob  ( 4 ,  0 ,  1 )  returns  exact  symbolic  probabilities,  and 
for  =9/10, 


^>(4,1)  = 

5  2  +  4  +1 

865000 

(  +1)5 

2476099 

ft(4,2)  = 

(5  2  +  4  +1) 

778500 

(  +1)5 

2476099 

fi>(4,3)  = 

2  (3  +1) 

29970 

(  +1)4 

130321 

^(4,4)  = 

3 

729 

(  +1)3 

6859 

0.3493 

0.3144 

0.2300 

0.1063, 


where  =  /  .  It  is  easy  to  verify  that  for  any  >  0,  i=l  /f)(4,  i)  =  1,  as  required. 

Example  2.  For  the  queue  described  in  Example  1,  calculate  the  fourth  flight’s  sojourn 
time  distribution,  mean  sojourn  time  and  sojourn  time  variance. 
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The  APPL  statements 


X  :=  ExponentialRV ( lambda) ; 
Y  :=  ExponentialRV (mu)  ; 

T  :=  Queue (X,  Y,  4,  0,  1)  ; 
Mean  (T) ; 

Variance (T) ; 


calculate  the  desired  results.  The  hist  two  lines  define  the  interarrival  and  service  time 
distributions,  while  the  third  line  calculates  the  fourth  bight’s  sojourn  time  distribution. 
The  last  two  lines  are  self  explanatory.  The  resulting  PDF  is 


f4(t)  =  4e~  r  (30  2  +  30  3r  +  24  +24  2  t  +  6  2  +  6  2  t  +  9t2  4+ 

12f2  3  +3 f2  2  2  +  t3  5 +  2 13  4  +t3  3  2)/6(  +  )5  t>  0. 

Using  f4{t)  above,  the  Mean  and  Variance  commands  return 

r  .  5  +  6  4  +  26  2  3  +  16  3  2  +  17  4  +  4  5 

E  ^  ; - + - 

(  +  ) 

and 

V  [74]  =  (181  2  8  +  484  3  7  +  816  4  6  +  868  5  5  +574  6  4+ 

244  7  3  +  40  9  +  68  8  2  + 12  9  +  10  +  4  10)  /  f  2(  +  )10)  . 

Substituting  =  1  and  =  10/9,  the  results  simplify  to 

Mi)  =  ^^e-10/9f(361r3+2109t2  +  5190t  +  5190)  t  >  0, 


E[Ta] 


23323347 

12380495 


1.8839,  and  V[T4] 


383506725720906  „ 

- ~  2.5021. 

153276656445025 


The  CPU  time  associated  with  the  examples  is  negligible.  Examples  1  and  2  represent  simple 
applications  of  these  procedures  that  circumvent  time  intensive  hand-calculations.  They  serve  only 
as  indications  of  more  challenging  problems  solvable  using  these  procedures. 

Example  3.  Calculate  the  mean  sojourn  time  of  the  30th  bight  in  an  M/M / 2  queue 
with  arrival  rate  =  1,  service  rate  =  9/20  (  =  10/9),  and  k  =  3  bights  initially 
present. 

The  mean  can  be  calculated  in  a  single  APPL  statement  by  embedding  the  function 
calls 

Mean (Queue (ExponentialRV ( 1 )  ,  ExponentialRV (9/20),  30,  3,  2)); 

which  yields  9.634524585. 


Representing  the  sojourn  time  distribution  for  the  nth  bight  in  closed  form  also  provides  valuable 
information  on  asymptotic  behavior  for  queueing  systems,  including  steady  state  convergence  rates 
for  different  initial  conditions.  Figure  1  shows  the  mean  sojourn  time  for  bight  n  =  1,2, . . . ,  120  in 
an  M/M/ 1  queue  with  =1,  =10/9,  and  =9/10  for  several  values  of  k.  The  points  that  are 

plotted  have  been  connected  by  lines.  As  expected,  despite  the  initial  condition,  all  cases  appear  to 
move  toward  the  steady-state  value  of  9  with  increasing  n.  The  horizontal  axis  is  only  limited  to 
n  =  120  for  display  purposes  and  in  fact,  identical  computations  were  earned  out  for  n  >  300  bights 
to  verify  convergence.  Flowever,  as  shown  in  the  cases  where  k  =  6  and  k  =  10,  the  convergence  to 
steady-state  is  not  always  monotone.  Additionally,  in  testing  various  traffic  intensities,  the  rate  of 
convergence  to  steady-state  increases  rapidly  with  decreasing  traffic  intensity  for  varying  values  of  k. 

APPL  also  has  the  ability  to  calculate  the  closed-form  cumulative  distribution  function  (CDF) 
for  the  nth  bight’s  sojourn  time  permitting  CDF  comparisons  for  varying  n  as  well  as  distribution 
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E[T„] 


Figure  1:  M/M/ 1  mean  sojourn  time  for  =  9/10  given  k  at  t  =  0. 


percentiles  for  a  given  flight.  Varying  k  for  an  M/M/ 1  queue  provides  another  basis  for  comparison 
of  CDFs,  where  resulting  CDFs  are  plotted  across  k  allowing  direct  comparison  of  sojourn  time 
percentiles  for  flight  n. 

Given  the  complete  specification  of  the  sojourn  time  distribution,  one  can  use  APPL  to  calculate 
not  only  the  mean  but  also  the  2nd,  3rd,  and  4th  moments  for  flight  n.  This  is  especially  valuable 
for  steady-state  analysis.  It  is  common  in  simulation  to  verify  attainment  of  steady-state  behavior 
by  examining  the  mean  delay  or  mean  sojourn  time.  Though  some  literature  exists  on  estimating 
transient  mean  and  variance,  we  are  not  aware  of  any  literature  addressing  higher  moments.  Literature 
addressing  the  second  moment  seems  mostly  focused  on  variance  estimation  and  not  necessarily 
convergence.  Therefore,  even  when  the  first  moment  might  acceptably  approximate  the  steady  state 
value,  there  is  reason  for  further  analysis  of  higher  moments.  For  example,  Figure  2  displays  the  first 
four  moments  of  the  sojourn  time  for  flight  n  in  an  M/M/ 1  queue,  where  =1,  =2,  =  1/2, 

with  the  initial  condition  k  =  0,4, 8. 

The  vertical  dashed  lines  give  the  smallest  flight  number  for  which  all  three  of  the  transient  values 
are  within  1%  of  the  steady  state  value.  The  relatively  low  traffic  intensity  =1/2  was  selected 
purposely  to  allow  quick  convergence  and  easy  visual  inspection.  Even  with  this  somewhat  low  traffic 
intensity,  it  is  apparent  that  the  higher  moments  converge  more  slowly  than  the  lower  moments.  In 
other  scenarios  where  >  1  /2,  the  higher  moments  exhibit  an  even  slower  convergence.  Each  vertical 
dashed  line  in  Figure  2  was  triggered  by  the  k  =  8  curve,  suggesting  that  the  moments  are  more  sensitive 
to  a  heavily  pre-loaded  system.  For  the  cases  k  =  0,4,8,  the  flight  numbers  for  which  the  transient 
results  were  within  1%  of  the  steady  state  values  are  listed  in  Table  3.  To  verify  the  initial-condition 
effect  on  the  convergence  rate  of  the  first  four  moments,  k  was  increasingly  incremented  beyond  eight 
and  displayed  a  further  slowing  of  convergence. 

Table  3:  Smallest  flight  number  where  the  sojourn  time  transient  result  is  within  1%  of  steady  state  for  an  M/M/1 
queue  with  £  =  0,4,8  and  =  1/2. 


k  =  0 

k  =  4 

k=  8 

E[T] 

19 

21 

36 

VVar[T] 

27 

29 

46 

E 

((T-  )/  ? 

28 

29 

50 

E 

({T-  )/  )4 

34 

35 

56 

5  FLIGHT  COVARIANCE  AND  CORRELATION  IN  THE  M/M/1  QUEUE 

The  dependence  exhibited  in  sojourn  times  of  successive  flights  is  one  reason  for  the  difficulty  in 
calculating  interval  estimators  for  queue  measures  of  performance.  In  the  simplest  case,  consider  an 
empty  and  idle  M/M/1  queue  with  interarrival  and  service  rates  and  .  Our  desire  is  to  calculate 
the  covariance  between  the  sojourn  times  of  any  two  flights  for  the  first  three  flights  of  the  day. 
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E[r„] 


VVar[7Y] 


E[((r„-|i)/o)3] 


E[((r„-|x)/o)4] 


Figure  2:  First  four  moments  of  the  M/M/1  sojourn  time  for  flights  2  through  60  for  =  1/2  and  k  =  0,4,8.  The 
arrival  rate  is  =1  and  the  service  rate  is  =  2,  resulting  in  steady-state  values  for  the  four  measures  of  performance 
of  1, 1,2,9. 


When  considering  n  =  3  flights,  there  are  five  possible  ways  flights  can  arrive  and  be  serviced. 
In  general,  for  n  flights,  the  number  of  ways  that  flights  can  enter  and  leave  the  ATC  is  given  by  the 
nth  Catalan  number,  which  is 

(2  n) ! 

"  (n!)(n  +  l)! 

Figure  3  shows  the  five  possible  arrangements  for  n  =  3  flights  along  with  the  sojourn  times  7j ,  73,  and 
73  for  each,  with  the  arrival  and  completion  times  for  the  ith  flight  denoted  by  a,-  and  respectively. 
The  vertical  arrows  at  event  times  represent  service  completions  (pointing  up)  or  arrivals  (pointing 
down). 
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< —  T2  — ►  < —  T3  — ►  Case  1 


c 


3 


Case  2 


c 


3 


Case  3 


Case  4 


c 


3 


Case  5 


Figure  3:  Five  cases  for  n  =  3  flights’  sojourn  times  in  an  M/M/1  queue. 


Using  a  theorem  presented  in  Kaczynski  et  al.  (2010),  the  joint  PDFs  for  each  of  the  pairs  (7j ,  73), 
(7) .  Tj),  and  ( 73 ,  Tj )  in  each  of  the  five  cases  can  be  determined  and  then  mixed  to  achieve  the  three 
associated  joint  PDFs.  Using  these  joint  densities,  the  symmetric  3x3  variance-covariance  matrix  is 

1  (2  +  )  2(  2  +  4  +5  2) 

2  (  +  )2  2  (  +  )4  2 

2  2  +4  +  2  (2  2  +  8  2  +  11  2  +2  3) 

(  +  )2  2  (  +  )4  2 

3  6  +  18  5  +45  4  2  +  54  3  3  _|_  30  2  4  +  g  5  +  6 


Substituting  =  1  and 

=  2,  for  example,  results  in 

'15  29  ' 

4  36  324 

7  13 

•  —  — 

'  0.2500  0.1389  0.0895  ' 
•  0.3889  0.2407 

18  54 

1451 

•  •  - 

L  2916  J 

•  •  0.4976 

These  results  have  been  verified  via  Monte  Carlo  for  the  first  n  =  3  flights.  The  sojourn  time 
variance  increases  with  flight  number  down  the  diagonal  of  the  matrix  because  of  the  nature  of  the 
queueing  process,  where  the  sojourn  time  distribution  for  each  additional  flight  is  dependent  on  all 
the  previous  flights.  On  the  other  hand,  the  off-diagonal  covariance  entries  in  each  row  decrease  with 
flight  separation,  for  example  [3  <  12. 
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6  EXTENDING  FLIGHT  COVARIANCE  CALCULATIONS 

The  APPL  procedure  Cov  (a,  b,  n)  calculates  the  covariance  between  flights  a  and  b  (a  <  b)  in 
a  system  of  n  flights.  For  computational  considerations  (i.e.,  evaluating  the  fewest  cases  necessary 
for  a  given  n),  setting  the  number  of  flights  n  =  b  provides  the  fastest  result.  Additionally,  calling 
Cov  ( a ,  b ,  n )  where  n  >  b  produces  a  result  identical  to  n  =  b  because  flights  arriving  after  flight  b  do 
not  affect  the  covariance  of  previous  flights.  Using  this  procedure,  the  symmetric  variance-covariance 
matrix  for  n  =  10  flights  with  parameters  =  1,  =  2,  and  =  1/2  is  calculated  and  displayed 

with  matrix  elements  rounded  to  four  significant  digits.  CPU  time  is  a  factor  in  these  computations. 
Each  element  in  the  tenth  column  of  the  variance-covariance  matrix  is  calculated  from  a  joint  PDF 
which  is  a  mixture  of  Cio  =  20! / (10!  11!)  =  16,796  component  distributions,  each  corresponding  to 
a  unique  ordering  of  ATC  arrivals  and  completions. 


'  0.2500 

0.1389 

0.0895 

0.0621 

0.0450 

0.0336 

0.0257 

0.0199 

0.0157 

0.0125 

• 

0.3889 

0.2407 

0.1639 

0.1176 

0.0872 

0.0663 

0.0513 

0.0402 

0.0319 

• 

• 

0.4976 

0.3251 

0.2286 

0.1676 

0.1263 

0.0972 

0.0759 

0.0600 

• 

• 

• 

0.5845 

0.3948 

0.2837 

0.2113 

0.1611 

0.1251 

0.0984 

• 

• 

• 

• 

0.6547 

0.4524 

0.3302 

0.2488 

0.1915 

0.1498 

• 

• 

• 

• 

• 

0.7119 

0.5000 

0.3694 

0.2808 

0.2177 

• 

• 

• 

• 

• 

• 

0.7587 

0.5396 

0.4022 

0.3080 

• 

• 

• 

• 

• 

• 

■ 

0.7974 

0.5725 

0.4298 

• 

• 

• 

• 

• 

• 

• 

• 

0.8293 

0.5999 

• 

• 

• 

• 

• 

• 

• 

• 

• 

0.8559 

As  the  traffic  intensity  increases,  so  do  the  values  in  the  variance-covariance  matrix.  To  illustrate, 
the  same  matrix  is  provided  for  the  increased  traffic  intensity  parameters  =1,  =  10/9,  and 

=  9/10.  The  increasing  sojourn-time  variance  along  the  diagonal  is  expected  with  the  increasing 
traffic  intensity.  In  addition,  the  rate  that  covariance  between  flights  decreases  as  flight  separation 
increases  is  less  pronounced. 


'  0.8100 

0.5856 

0.4737 

0.4040 

0.3553 

0.3189 

0.2904 

0.2673 

0.2481 

0.2318 

• 

1.3956 

1.1097 

0.9393 

0.8226 

0.7363 

0.6692 

0.6150 

0.5702 

0.5323 

• 

• 

1.9561 

1.6298 

1.4167 

1.2626 

1.1441 

1.0494 

0.9714 

0.9057 

• 

• 

• 

2.5021 

2.1458 

1.8995 

1.7142 

1.5679 

1.4484 

1.3485 

• 

• 

• 

• 

3.0364 

2.6565 

2.3831 

2.1715 

2.0009 

1.8593 

• 

• 

• 

• 

• 

3.5605 

3.1614 

2.8652 

2.6310 

2.4389 

• 

• 

• 

• 

• 

• 

4.0754 

3.6600 

3.3444 

3.0904 

• 

• 

• 

• 

• 

• 

• 

4.5818 

4.1524 

3.8199 

• 

• 

• 

• 

• 

• 

■ 

• 

5.0803 

4.6386 

• 

• 

• 

• 

• 

• 

• 

• 

• 

5.5713 

Traditional  steady-state  queueing  theory  and  analysis  lacks  the  insight  provided  in  these  transient 
variance-covariance  matrices.  For  military  airfields  where  the  number  of  flights  in  a  day  is  so  small  that 
true  steady  state  is  never  achieved,  routine  queueing  measures  of  performance  are  not  representative 
of  reality.  Additionally,  consider  a  system  where  the  traffic  intensity  exceeds  one.  An  appropriate 
example  is  the  Haiti  airport  following  the  2010  earthquake  where  military  control  over  the  airport 
was  exercised.  For  a  such  a  system,  an  analyst  might  be  interested  in  the  covariance  between  flight 
sojourn  times.  Increasing  the  traffic  intensity  so  that  >  1  does  not  preclude  covariance  calculations 
using  this  method,  and  therefore  allows  transient  analysis  of  such  systems.  A  variance-covariance 
matrix  for  =3/2,  is  presented  below.  Given  this  traffic  intensity,  the  system  is  unstable  and  the 
expected  sojourn  times  for  successive  flights  increase  without  bound.  Along  the  main  diagonal  the 
flight  variance  is  clearly  increasing,  and  the  covariance  decreases  as  the  separation  occurs  between 
flights.  This  decrease  is  monotonic,  and  though  not  studied  in  detail  here,  it  appears  that  the  rate  of 
covariance  decrease  might  be  of  interest  for  an  unstable  traffic  intensity. 
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.2500 

1.8900 

1.7172 

1.6135 

1.5438 

1.4937 

1.4558 

1.4263 

1.4027 

1.3835 

• 

4.1400 

3.7368 

3.5018 

3.3459 

3.2344 

3.1507 

3.0856 

3.0337 

2.9913 

• 

• 

6.0957 

5.6825 

5.4166 

5.2292 

5.0896 

4.9817 

4.8958 

4.8261 

• 

• 

• 

8.1312 

7.7208 

7.4397 

7.2332 

7.0747 

6.9493 

6.8479 

• 

• 

• 

• 

10.2424 

9.8410 

9.5538 

9.3361 

9.1652 

9.0276 

• 

• 

• 

• 

• 

12.4235 

12.0342 

1 1 .7463 

11.5230 

11.3444 

• 

• 

• 

• 

• 

• 

14.6687 

14.2931 

14.0081 

13.7828 

• 

• 

• 

• 

• 

• 

• 

16.9727 

16.6115 

16.3319 

• 

• 

• 

• 

• 

■ 

• 

• 

19.3310 

18.9846 

• 

• 

• 

• 

• 

• 

• 

• 

• 

21.7397 

7  SOJOURN  TIME  COVARIANCE  WITH  k  FLIGHTS  INITIALLY  PRESENT 


When  k  flights  are  present  in  the  M/M/ 1  queue  at  time  zero,  the  approach  used  to  compute  sojourn¬ 
time  covariance  between  flights  becomes  more  difficult.  When  the  two  flights  of  interest  possess 
indices  larger  than  k  (i.e.,  T,  where  i  >  k).  then  the  approach  is  similar  to  that  derived  in  Section  6  of 
Kaczynski  et  al.  (2010).  However,  there  are  two  other  possibilities.  The  first  possibility  is  that  the 
first  flight  has  an  index  of  k  or  less,  and  the  second  flight  has  an  index  larger  than  k.  In  this  instance, 
the  only  difference  in  deriving  the  joint  CDF  is  that  the  lower  indexed  flight  begins  its  sojourn  time 
at  time  zero.  In  the  second  possibility,  both  flights  have  an  index  of  k  or  below.  If  these  indices  are  i 
and  j,  where  i  <  j  <  k,  the  time  intervals  for  sojourn  times  T,  and  Tj  begin  at  zero.  These  calculations 
are  coded  in  Maple  as  the  procedure  kCov  (X,  Y,  a,  b,  n,  k) .  The  first  two  arguments  X  and 
Y  are  the  distribution  of  time  between  arrivals,  exponential  ),  and  the  service  time  distribution, 
exponential  ),  respectively.  The  arguments  a  and  b  are  the  flights  of  interest  for  the  covariance 
calculation,  where  a  <b.  The  argument  n  is  the  number  of  flights  processing  through  the  system  not 
including  those  present  at  time  zero,  which  is  indicated  by  the  last  argument,  k.  Therefore,  the  total 
number  of  flights  processing  through  the  system  is  n  +  k,  and  a  covariance  calculation  between  any 
two  of  these  flights  can  be  achieved  with  the  appropriate  function  call.  For  example,  the  function  call 
kCov  (ExponentialRV  ( 1 )  ,  ExponentialRV  ( 2 )  ,  1,  2,  1,  3)  calculates  the  covariance 
between  flights  1  and  2  in  an  M/M/1  queue  with  an  arrival  rate  =  1,  with  service  time  rate  =  2, 
with  three  flights  present  at  time  zero,  and  a  single  additional  flight  processes  through  the  system. 
The  variance-covariance  matrix  for  an  M/M/1  queue  with  arrival  rate  =  1,  and  service  time  rate 
=  2,  where  k  =  4  flights  are  present  at  time  zero  and  an  additional  n  =  6  flights  process  through 
the  system  is 


4 


263490131  208262483  4506205633 


172186884 


172186884 

63939878 

43046721 


4649045868 

1359189250 

1162261467 

179260456277 

125524238436 


l 

1 

211 

1579 

11651 

28553 

630131 

4646155 

4 

4 

972 

8748 

78732 

236196 

6377292 

57395628 

1 

1 

211 

1579 

11651 

28553 

630131 

4646155 

2 

2 

486 

4374 

39366 

118098 

3188646 

28697814 

3 

3 

211 

1579 

11651 

28553 

630131 

4646155 

4 

4 

324 

2916 

26244 

78732 

2125764 

19131876 

1 

211 

1579 

11651 

28553 

630131 

4646155 

243 

2187 

19683 

59049 

1594323 

14348907 

37289 

271153 

1966777 

1588153 

34755203 

763875281 

26244 

236196 

2125764 

2125764 

57395628 

1549681956 

1629655 

11663887 

9353743 

203800469 

4465399991 

1062882 

9565938 

9565938 

258280326 

6973568802 

98323535707 

125524238436 

29402061622 

31381059609 

379721786263 

3389154437772 

62708955663745 

45753584909922 


Here  the  matrix  elements  are  listed  exactly,  illustrating  Maple/ APPL’s  ability  to  conduct  exact  cal¬ 
culations.  Unlike  the  previous  variance-covariance  matrices,  some  row  elements,  in  particular  those 
elements  associated  with  flights  that  are  initially  present,  do  not  decrease  monotonically.  These  entries 
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are  explained  in  detail  in  Kaczynski  et  al.  (2010).  It  can  be  shown  in  general  that 

Var(2/)=Cov(2/,rJ-)  =  -T 


for  i  <  j  <  k,  where  k  flights  are  present  at  time  zero. 

8  OTHER  POTENTIAL  APPLICATIONS 

Queueing  applications  in  which  the  queue  behavior  does  not  reach  steady-state  are  commonplace 
in  the  military.  One  potential  area  of  study  is  the  soldier  processing  center  in  Kuwait.  This  center 
is  a  clearing  house  for  processing  service  members  in  and  out  of  Iraq  and  Afghanistan.  The  entire 
process  is  complicated,  however,  individual  stations  certainly  follow  the  FIFO  queue  discipline  and 
are  good  candidates  for  Markovian  queues.  This  center  is  known  to  have  substantial  problems  in 
(1)  arrival  rates  exceeding  service  rates,  i.e.,  >  1,  (2)  stations  that  shut  down  intermittently  without 

servicing  all  customers  already  present,  i.e,  do  not  achieve  steady  state,  and  (3)  form  queues  prior 
to  servers  being  available  for  service,  i.e.,  k  >  0  customers  present  at  time  t  =  0.  Another  potential 
application  is  installation  access  through  security  checkpoints.  It  is  often  observed  that  the  arrival  rate 
is  substantially  higher  than  the  service  rate,  especially  during  peak  periods.  This  is  better  represented 
as  a  non-homogeneous  Poisson  process,  and  would  be  an  interesting  future  research  topic  in  transient 
queueing  analysis. 

9  CONCLUSIONS  AND  LIMITATIONS 

Previous  transient  analysis  results  for  the  M/M/1  and  M/M /s  queues  have  been  combined  with  the 
functionality  of  the  Maple  computational  engine  (and  subsequently  APPL)  in  a  military  air  traffic 
control  environment  to  develop  both  symbolic  and  numeric  exact  flight  sojourn  time  PDFs  that  can  be 
manipulated  to  compute  and  study  various  measures  of  performance.  A  complete  variance-covariance 
matrix  for  the  first  n  =  10  flights  and  varying  traffic  intensity  is  calculated,  illustrating  this  approach’s 
ability  to  determine  the  joint  PDF  between  two  flight  sojourn  times.  The  models  presented  in  this 
paper  have  assumed  that  there  are  two  runways  that  are  dedicated  to  incoming  and  outgoing  flights 
which  operate  independently.  Therefore,  both  arriving  flights  and  departing  flights  can  be  viewed 
separately  as  independent  M/M/1  queues,  with  the  runway  playing  the  role  of  the  server.  An  analyst 
would  need  to  assure  that  data  supports  the  Markovian  assumptions,  and  the  policy  at  the  airport 
supports  the  FIFO  queue  discipline.  A  single-runway  airport  handling  both  arriving  and  departing 
flights  is  less  likely  to  conform  to  the  models  presented  here.  A  logical  epoch  for  arrival  time  for 
an  outgoing  flight  would  be  the  hand  off  from  ground  control  to  tower  for  outgoing  flights;  a  logical 
epoch  for  arrival  time  for  an  incoming  flight  would  be  initial  report  to  terminal  control.  Similar 
conventions  would  apply  when  determining  the  appropriate  time  for  completion  of  service.  The 
results  offer  a  framework  for  describing  how  the  well-known  M/M /s  queue  steady-state  results  occur 
as  the  queue  progresses  toward  steady-state.  When  possible,  results  are  checked  against  corresponding 
Monte  Carlo  simulation  and/or  previous  literature.  The  first  principles  derivation  suggests  a  viable 
alternative  for  future  research  would  be  to  apply  the  approaches  provided  in  this  work  to  a  G/G/l 
queue.  Computational  considerations  take  priority  as  n  increases.  Making  use  of  other  computational 
formulae  (such  as  Hagwood  (2009))  may  offer  significant  time  savings  and  is  another  interesting 
avenue  for  future  work. 
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